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Abstract 

In  this  paper,  an  approximate  technique  is  used  to  obtain  pressure  and  velocity  distribution  in  the  cathode  chamber  of  PEM  fuel  cells.  The 
technique  is  based  on  polynomial  profile  approximation.  The  technique  predicts  two  dimensional  pressure  and  velocity  distribution  in  closed 
form  as  a  function  of  independent  variables.  The  approximate  solution  developed  compares  reasonably  with  the  rigorous  numerical  solution. 
©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

Fuel  cells  are  devices  that  convert  chemical  energy  di¬ 
rectly  into  electrical  energy.  Fuel  cells  are  believed  to  be  a 
most  promising  power  source  for  a  wide  range  of  applications 
by  virtue  of  their  high-energy  efficiency,  pollution  free  char¬ 
acteristics,  compactness  in  design  and  operation.  In  the  last 
decade  significant  improvement  has  been  made  in  the  design 
of  PEM  fuel  cells. 

Models  for  PEM  fuel  cells  have  features  both  com¬ 
mon  and  different  from  models  for  secondary  batteries. 
Because  of  recent  interest  in  PEM  fuel  cells,  researchers 
from  various  disciplines  have  started  using  different  kind 
of  models  for  PEM  fuel  cells  according  to  their  conve¬ 
nience/objectives/understanding  [1-15].  Recently  new  mod¬ 
els  based  on  thermodynamics  related  equations  have  been 
developed  to  predict  the  performance  of  PEM  fuel  cells  [28]. 
In  order  to  accurately  predict  the  performance  of  PEM  fuel 
cells  one  has  to  solve  different  dependent  variables  including 
pressure,  concentration  (or  mole  fraction),  electrolyte  poten¬ 
tial,  temperature  etc  [29,30].  Models  varying  from  1D-3D, 
one  region  to  3  or  5  regions  have  been  used  in  the  litera¬ 
ture  [16-31].  Even  for  predicting  two-dimensional  pressure 
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or  velocity  distribution  accurately,  only  numerical  solution 
has  been  used  and  reported  in  the  literature. 

Most  of  the  model  equations  for  PEM  fuel  cells  are  non¬ 
linear.  However,  the  pressure  distribution  for  oxygen  inside 
the  cathode  chamber  (Fig.  1)  can  be  sufficiently  described  us¬ 
ing  the  Laplace’s  equation  [25].  Though  the  governing  equa¬ 
tion  is  linear  (Laplace’s  equation),  researchers  haven’t  used 
analytical  solution  for  pressure  as  the  boundary  condition  at 
the  top  ( y  =  H  in  Fig.  1)  has  three  different  boundary  con¬ 
ditions  depending  on  the  value  for  x  [25,26].  In  this  paper, 
we  arrive  at  approximate  solutions  for  pressure  and  velocity 
distributions  [32,34,35].  These  approximate  solutions  will 
be  fed  into  the  non-linear  convective-diffusion  equation  for 
concentration  distribution  in  a  later  publication  to  arrive  at 
analytical/efficient  solution  for  two-dimensional  concentra¬ 
tion  profiles  for  oxygen  inside  the  cathode  chamber  shown 
in  Fig.  1.  The  approximate  model  developed  is  compared 
with  a  rigorous  numerical  solution  to  analyze  the  accuracy  of 
approximation. 

1.1.  Approximate  models  for  pressure  and  velocity 
distributions 

The  cathode  chamber  shown  in  the  Fig.  1  is  modeled  in 
this  paper.  The  following  assumptions  are  made: 
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Fig.  1 .  Schematic  of  the  cathode  chamber  modeled. 


•  Water  exists  only  in  the  vapor  form  and  the  effects  of  pore 
plugging  are  neglected  [42,43]. 

•  The  fuel  cell  is  isothermal. 

•  Darcy’s  law  is  valid. 

•  Fluid  is  incompressible. 

•  Temperature  does  not  affect  pressure. 

•  The  electrode  layer  is  assumed  to  be  homogeneous  with 
uniform  physical  properties. 

•  Transport  of  the  reactants  in  the  electrode  layer  is  ne¬ 
glected. 


The  governing  equation  for  the  pressure  is  given  by 
Laplace’s  Eq.  [25] 


73^.y))  +  l^p(*.y))  =  o 


dx 2 


dy2 


(1) 


The  boundary  conditions  are  summarized  as  follows. 


At  A  =  0, 
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At  A  =  L, 
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At  the  membrane  (y  =  0), 

dp 

dy 

=  0  (0  <  x  <  L) 

(4) 

At  the  inlet, 

P  =  Pin 

(y  = 

H  and  0  <  x  <  L  i ) 

(5) 

At  the  current  collector. 


dp 

dy 


=  0 


(y  =  H  and  L\  <  x  <  Lj) 


(6) 


At  the  outlet,  p  —  pout  (y  =  H  and  Lj  <  x  <  L)  (7) 

where,  pm  is  the  inlet  pressure  and /?out  is  the  outlet  pressure. 
The  following  dimensionless  variables  are  introduced, 


where  e  =  ///Lis  the  aspect  ratio.  The  dimensionless  boundary 
conditions  are  given  by. 


At  A  =  0, 


0  <  Y  <  1 


(10) 


AtA  =  1, 


0  <  Y  <  1 


(11) 


At  the  membrane  (T  =  0), 


(0  <  A  <  1) 

(12) 


At  the  inlet  P  —  1 


(V  =  1  and  0  <  A  < 


(13) 


At  the  current  collector 


L  i  Li\ 
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At  the  outlet  P  —  0 


L? 

1  and  —  <  A  <  1 
L  ~  ~ 


(15) 


Next,  a  polynomial  approximation  is  made  for  pressure  in 
each  region  [11,12].  For  region  1  pressure  is  taken  as, 

Pl(X,Y)  =  a0(X)  +  al(X)Y2  (16) 


Forregion2,  Pj{ A,  Y)  =  bo(X )  +  b[(X)Y2  (17) 

and  for  region  3,  P^(  A,  Y)  =  co(A)  +  ci(A)T2  (18) 

Note  that  the  term  ‘F  is  left  in  Eqs.  ( 16)— ( 18)  to  sat¬ 
isfy  the  boundary  condition  at  Y=0  (Eq.  (12)).  Next,  the 
functions  «o(A),  ai(A),  bo(X),  b[(X ),  co(A)  and  ci(A)  are 
found. The  volume  averaged  equation  for  pressure  in  the 
first  region  (7J|  ave(A))  is  given  by  integrating  the  approx¬ 
imation  equation  for  P\(X,Y)  in  y  direction  from  0  to 
1. 

r 1  i 

/  Pi  (A,  Y)  dv  =  Piave(A)  =  a0(A)  +  -aj(A)  (19) 
Jo  4 

By  applying  the  boundary  condition  at  Y=  1  in  the  first  region 
i.e.,  substituting  Eq.  (13)  in  Eq.  (16)  we  get: 


P  = 


P  ~  Lout 
Pin  Pout 


y=A 

H 


x 

A  =  - 
L 


(8) 


The  governing  equation  and  boundary  conditions  in  dimen¬ 
sionless  form  are. 


d1 

dX2 


P(  A,  T))  +  (^P(A,  T)  )  =  0 


(9) 


a0(A)  +  ai(A)  —  1  =0 

Eqs.  (17)  and  (18)  are  solved  to  get 

a0(A)=  ^Piave(A)-  i 

flt(A)  =  — ^Pjave(A)  +  ^ 


(20) 


(21) 
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By  substituting  Eq.  (21)  in  Eq.  (16)  we  get, 

MX,  7)  =  ^Piave(X)  -  l-  +  ^Piave(X)  +  0  72 

(22) 

Next,  the  approximation  equation  (Eq.  (22))  is  substituted  in 
the  governing  equation  (Eq.  (8))  and  integrated  in  Y  from  0 
to  1  to  get 

—  3P]'dve(X)  +  3  =  0  (23) 

Similarly  the  governing  equations  for  volume  averaged  pres¬ 
sures  in  regions  2  and  3  are  obtained  after  similar  steps 
as: 

g2  P2ave(z)) = °  (24) 


2  fT^piave(Z)^) 


\dX2 


s 


2 


P3  ave(X) 


-  3  P3 ave(X)  =  0 


(25) 


Note  that  a  single  partial  differential  equation  (Eq.  (9))  is  con¬ 
verted  to  three  second  order  ordinary  differential  equations 
(Eqs.  (23),  (24),  and  (25)).  We  need  six  boundary  conditions 
to  solve  this  system  of  boundary  value  problems.  We  get  two 
boundary  conditions  by  volume  averaging  the  boundary  con¬ 
ditions  at  X=  0  and  X=  1  (Eqs.  (10)  and  (11)). 


d  P\  ave 
dX 


=  0  at  X  =  0 


(26) 


d  P3  ave 
dX 


=  0  at  X  =  1 


(27) 


Next,  the  pressure  and  its  gradient  should  be  continuous  at 
X=  ( 1/4)  and  at  7  =  (3/4).  By  volume  averaging  the  continuity 
of  pressure  and  its  gradient  we  obtain, 

AtX  =  Piave(X)  =  P2  ave(X)  (28) 


dPiave(X)  dP2ave(X) 
dX  ~~  dX 


(29) 


and  at  X  =  |,  P2ave(X)  =  P3ave(X)  (30) 

dP2ave(X)  _  dP3ave(X) 
dX  ~  dX 

Solving  the  governing  equations  (Eqs.  (23),  (24)  and  (25)) 
using  the  six  boundary  conditions  given  in  the  Eqs.  (26), 
(27),  (28),  (29),  (30)  and  (31)  we  get, 


Piave(X)  =  1  — 


2s  cos  h(\/3X/s) 

4  cos  h(^/3/4s)  +  sin  h(\/3/4s)V3 


(32) 


P2  ave(X)  =  1  + 


2  sin  h(V3/4s)V3X 
4  cos  h(4 3/4e)  +  sin  h(y/3/4s)V3 


1  sin/j(V3/4e)V3 

2  4  cos  h(y/ 3 /4s)  +  sin  h(j3/4s)V3 


2  cos  /z(V3/4e) 

4  cos  h(\/3/4s)  +  sin  /i(V3/4e)V3 


(33) 


P3ave(X)  = 


2s  cos/?(\/3(l  —  X)/s) 

4  cos  /?(V3/4e)  +  sin  /;(V3/4e)  V3 


(34) 


By  substituting  the  above  expressions  in  Eqs.  (16),  (17),  (18) 
and  (22)  we  get 


Pi(X,  F)=  1- 


3s  cos  h(y/3X/s) 

4  cos  h(V3/4s)  +  sin  h(V3/4s)V3 


3s  cos  li(y/3X/s) 

4  cos  h(y/3/4s)  +  sin  h(^/3/4s)j3 


Y2 

(35) 


MX,  Y)=  1  - 


2  sin  h(V3/4s)V3 
4  cos  h(\Z3/4s)  +  sin  h(\/3/4s)V3 


1  sin  h(V3/4s)*j3 

2  4  cos  h(y/3/4e)  +  sin  h(V3/4s)V3 

2s  cos  h(V3/4s) 

4  cos  h(V3/4s)  +  sin  h(V3/4s)V3 


MX,  Y)  = 


3s  cos/i(\/3(l  —  X)/s) 

4  cos  h(\/3/4s)  +  sin  /i(V3/4e)V3 


3scosh(V3X/s)  y2 

4  cos  h(j3/4s)  +  sin/?(V3/4e)V3 

(37) 


Note  that  a  closed  form  solution  for  pressure  is  obtained  in 
terms  of  aspect  ratio  s.  The  dimensionless  velocity  compo¬ 
nents  U  and  V in  x  and  y  directions  respectively  are  calculated 
using  the  following  equations 


U  = 

V  = 


9P 

'dX 

9P 

dY 


(38) 

(39) 


2.  Results  and  discussion 

Parameters  (L\,  L2,  L  and  e)  in  the  model  are  taken  from 
literature  [25]  and  are  listed  in  Table  1.  Fig.  2  shows  the 
distribution  for  dimensionless  pressure  predicted  using  the 
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Table  1 

Modeling  parameters  in  the  electrode 

Cathode  width,  L 

0.18  cm 

Cathode  height,  H 

0.02  cm 

Length  of  region  1,  L\ 

L/4 

Length  of  region  2,  L2 

LI  2 

Length  of  region  3,  L3 

L/4 

approximate  model  developed.  From  Fig.  2  we  observe  that 
pressure  varies  linearly  in  the  cathode  chamber  between  the 
current  collector  and  the  membrane. 

In  the  following  figures,  rigorous  numerical  solution  ob¬ 
tained  by  solving  Eq.  (9)  with  Eqs.  (10)-(  14)  using  FEMLAB 
[32]  is  used  for  comparison. 

Fig.  3,  shows  the  profiles  for  dimensionless  pressure  at 
the  boundary  X=  0.  The  approximate  solution  developed  in 
this  paper  is  compared  with  a  rigorous  2D  numerical  solution 


x  y 

Fig.  2.  Profile  for  dimensionless  pressure  distribution  in  the  cathode  of  a 
PEM  Fuel  cell. 


Fig.  3.  Dimensionless  pressure  distribution  at  X = 0.  Dotted  line  corresponds 
to  the  approximate  model  developed.  Solid  line  corresponds  to  the  rigorous 
numerical  solution. 


Y 


Fig.  4.  Dimensionless  pressure  distribution  atX=  1 .  Dotted  line  corresponds 
to  the  approximate  model  developed.  Solid  line  corresponds  to  the  rigorous 
numerical  solution. 

obtained  using  FEMLAB,  a  finite  element  software  [32] .  Both 
the  curves  started  at  0.993  and  ended  at  1. 

Fig.  4  shows  the  profile  for  dimensionless  pressure  at  the 
boundary  X  =  1 .  In  this  case  also  there  is  a  minor  difference 
in  the  magnitude  (around  0.001). 

Fig.  5  shows  the  profile  for  dimensionless  pressure  at  the 
boundary  Y=  0.  At  this  boundary  cathode  channel  and  mem¬ 
brane  come  into  contact.  Knowledge  of  gas  pressure  at  this 
interface  is  important  because,  oxygen  passes  through  cath¬ 
ode  gas  channel  to  the  cathode  catalyst  layer  where  the  elec¬ 
trochemical  reaction  takes  place.  The  curves  using  numerical 
solution  and  the  approximate  model  developed  matches  accu¬ 
rately  excepting  that  there  is  slight  bump  in  the  curve  obtained 
using  maple  at  X=  3/4,  this  could  be  because  of  the  reason 
that  at  the  interface  between  second  and  third  regions  i.e.,  at 
X  =  3/4,  average  pressures  and  their  differentials  are  equated 
instead  of  equating  pressures  at  values  of  X. 


Fig.  5.  Dimensionless  pressure  distribution  at  the  membrane  interface 
(F=0).  Dotted  line  corresponds  to  the  approximate  model  developed.  Solid 
line  corresponds  to  the  rigorous  numerical  solution. 
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Fig.  6.  Dimensionless  pressure  distribution  at  Y=  1 .  Dotted  line  corresponds 
to  the  approximate  model  developed.  Solid  line  corresponds  to  the  rigorous 
numerical  solution. 


Fig.  6  shows  the  profile  for  dimensionless  pressure  at  the 
boundary  Y=  1 .  Both  the  curves  using  FEMLAB  and  the  ap¬ 
proximate  model  developed  matches  well  in  this  case,  ex¬ 
cepting  at  the  interface  between  first  and  second  regions  i.e., 
at  X=  1/4.  The  reason  for  this  anomaly  could  be  same  as 
mentioned  in  the  Fig.  5. 

Figs.  7  and  8  show  the  comparison  of  vertical  components 
of  dimensionless  velocity  (V).  The  approximate  model  de¬ 
veloped  predicts  the  behavior  reasonably  well. 

Fig.  9  shows  the  comparison  of  profiles  for  X  component 
of  dimensionless  velocity  U  at  Y=  0.  The  approximate  model 
developed  fails  in  capturing  the  behavior,  especially  at  the  in¬ 
terfaces  at  X  =  1/4  and  X  =  3/4.  To  improve  the  accuracy  one 


Y 


Fig.  7.  Comparison  of  y  component  of  dimensionless  velocity  (V)  distribu¬ 
tion  atX=0  and  Y  from  0  to  1.  Dotted  line  corresponds  to  the  approximate 
model  developed.  Solid  line  corresponds  to  the  rigorous  numerical  solution. 


Fig.  8.  Comparison  of  y  component  of  dimensionless  velocity  (V)  distribu¬ 
tion  at  X=  1  and  Y  from  0  to  1.  Dotted  line  corresponds  to  the  approximate 
model  developed.  Solid  line  corresponds  to  the  rigorous  numerical  solution. 


=> 


Fig.  9.  Comparison  of  x  component  of  dimensionless  velocity  ( U)  distribu¬ 
tion  at  Y  =  0.  Dotted  line  corresponds  to  the  approximate  model  developed. 
Solid  line  corresponds  to  the  rigorous  numerical  solution. 


=) 


X 

Fig.  10.  Comparison  of  x  component  of  dimensionless  velocity  ( U )  distribu¬ 
tion  at  Y  =  0.  Dotted  line  corresponds  to  the  approximate  model  developed. 
Solid  line  corresponds  to  the  rigorous  numerical  solution.  Dashed  line  cor¬ 
responds  to  the  improved  approximate  model. 
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can  add  more  number  of  terms  in  the  polynomial  approxi¬ 
mation  (Eqs.  (16)  to  (18)).  Alternately  we  found  that  if  we 
define 

Piave(X)  =  f  f(y)P]  (x,  y)dy 
Jo 

and  redo  the  calculations,  we  obtain  better  results  as  shown 
in  Fig.  10.  We  are  currently  working  on  optimizing  a  function 
f(y )  which  can  be  used  to  obtain 

Pi ave(X)=f  f{y)P\(x,  y)dy 

Jo 

we  plan  to  publish  this  later. 

3.  Conclusions 

In  this  paper,  an  approximate  technique  was  used  to  obtain 
closed  form  solution  for  pressure  and  velocity  distributions 
[32-41]  as  a  function  of  system  parameter,  the  aspect  ratio 
e.  The  approximate  model  predicts  the  behavior  reasonably 

well. 

For  modeling  PEM  fuel  cells  in  hybrid  environment  one 
has  to  simulate  models  for  PEM  fuel  cells  in  series/parallel 
combination  with  models  for  batteries  and  other  electrical 
components.  This  paper  is  our  first  step  towards  simplifying 
PEM  fuel  cell  models  for  hybrid  applications.  In  our  future 
publications,  we  plan  to  use  the  approximate  model  devel¬ 
oped  for  pressure  and  velocity  distribution  to  solve  the  con¬ 
vection  -  diffusion  equation  [19,20,40,41]  for  oxygen  con¬ 
centration  distribution  in  closed-form  efficiently. 
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